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ABSTRACT 

To critically assess the binary compact object merger model for short gamma ray bursts (GRBs) - 
specifically, to test whether the short GRB rates, redshift distribution and host galaxies are consis- 
tent with current theoretical predictions - it is necessary to examine models that account for the 
high-redshift, heterogeneous universe (accounting for both spirals and ellipticals). We present an in- 
vestigation of predictions produced from a very large database of first-principle population synthesis 
calculations for binary compact mergers with neutron stars (NS) and black holes (BH), that sample 
a seven-dimensional space for binaries and their evolution. We further link these predictions to (i) 
the star formation history of the universe, (ii) a heterogeneous population of star-forming galaxies, 
including spirals and ellipticals, and (iii) a simple selection model for bursts based on flux-limited 
detection. We impose a number of constraints on the model predictions at different quantitative 
levels: short GRB rates and redshift measurements, and, for NS-NS, the current empirical estimates 
of Galactic merger rates derived from the observed sample of close binary pulsars. Because of the 
relative weakness of these observational constraints (due to small samples and measurement uncer- 
tainties), we find a small, but still substantial, fraction of models are agreement with available short 
GRB and binary pulsar observations, both when we assume short GRB mergers are associated with 
NS-NS mergers and when we assume they are associated with BH-NS mergers. Notably, we do not 
need to introduce artificial models with exclusively long delay times. Most commonly models produce 
mergers preferentially in spiral galaxies, in fact predominantly so, if short GRBs arise from NS-NS 
mergers alone. On the other hand, typically BH-NS mergers can also occur in elliptical galaxies (for 
some models, even preferentially), in agreement with existing observations. As one would expect, 
model universes where present-day BH-NS binary mergers occur preferentially in elliptical galaxies 
necessarily include a significant fraction of binaries with long delay times between birth and merger 
(often O(lOGyr)); unlike previous attempts to fit observations, these long delay times arise naturally 
as properties of our model populations. Though long delays occur, almost all of our models (both a 
priori and constrained) predict that a higher proportion of short GRBs should occur at moderate to 
high redshift (e.g., z > 1) than has presently been observed, in agreement with recent observations 
which suggest a strong selection bias towards successful follow-up of low-redshift short GRBs. Finally, 
if we adopt plausible priors on the fraction of BH- NS mergers with approp riate combination of spins 
and masses to produce a short GRB event based on lBelczvnski et ail (|2007f ). then at best only a small 
fraction of BH-NS models could be consistent with all current available data, whereas NS-NS models 
do so more naturally. 

Subject headings: Stars: Binaries: Close; Gamma-ray bursts 



1. INTRODUCTION 

The afterglows of several short-hard gamma ray bursts 
(GRBs) have recently been localized on the sky, al- 
lowing reasonably precise determina tion of their hosts , 
redshif t s, an d energetics (see, e.g.. i Berger et all 120061 : 
iBergerl 120061 : iBarthelmv et all 12003 : iFox et all l2005f h 
The energetics, presence in both old and star form- 
ing host galaxies, absence of supernovae afterglow char- 
acteristics, and in some cases apparent host offsets 
of these bursts seem qualitatively consistent with the 
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merger hypothesis (MH): the notion that most short- 
hard bursts arise from the disruption of a neutron 
star (NS) in either a NS -NS or black hole-NS binary 
(see, e.g, iLee et all (|2005| ) and refe rences ther e in; th is 
GRB model was first presented by iPaczvnskl (|1986h ). 
While other populations, such as bursts from magne- 
tars, may contribute to the total short GRB event 
rate, the fraction of GRBs produced by nearby mag- 
netars is believed to be small (see, e.g.. | Nakar etaLI 
l2006tlPopov & Sternll2006l : lLazzati et al.ll2005f) . Further- 
more, empirical and theoretical estimates for compact 
object merger rates based on studies of the Milky Way 



(see, e . g.. iKim et ai.i i^uuat lU'snaugnnessv et ai.i izuua. 
2007ri iBelczvnski et all I2002U l2006bl: iKaloeera et all 



e. g.. iKim et all 12001 10'ShaughnessvlTaTl I2005L 



2004; fNagamine et al.ll2006f ) are r oughly consistent with 



BATSE and Swi ft observations (jGuetta fc Piranl 120061 . 
[2lMlAlIdol[200l . 

The number of observed radio pulsars with neutron 
star companions can provide a robust quantitative test 
of the MH. For example, using well-understood selec- 
tion effects and fairly minimal population modeling (i.e., 
a luminosi t y func tion and a beaming correction factor), 
IKim et all (|2003f ) developed a statistical method to de- 
termine which double neutron star coalescence rates were 
consistent with NS-NS seen in the Milky Way. However, 
in distinct contrast to NS-NS in the Milky Way, little is 
known directly about the short GRB spatial or luminos- 
ity distribution. 

Short GRBs are still detected quite infrequently (i.e, a 
handful of detections per year for Swift); sufficient statis- 
tics are not available for a robust nonparametric estimate 
of their distribution in redshift z and peak luminosity 
L. To make good use of the observed (z, L) data, we 
must fold in fairly strong prior assumptions about the 
two-dimensional density d 3 N/dtdLdz(L,z). Typically, 
these priors are constructed by convolving the star for- 
mation history of the universe with a hypothesized dis- 
tribution for the "delay time" between the short GRB 
progenitor's birth and the GRB event, as well as with 
an effective (detection- and angle-averaged) luminosity 
distribution for short GRBs. Observations are thus in- 
terpreted as constraints on the space of models, rather 
than as direct measurements of the z,L distribution 
(lAnddl2004l:lGuetta fc Piranl[2005l [20061 : 1 Gal- Yam et all 
|2005| ). A similar technique has been applied with con- 
siderable success to long GRB observations (see, e.g., 
Porciani fc Madaull2001l : iGuetta k Piranll2005t ISchmidtl 



19991 : IChe et al.l I1999I . and references therein) : as ex- 
pected from a supernovae origin, the long GRB rate is 
consistent with the star formation history of the uni- 
verse. And within the context of specific assumptions 
about the merger delay time distribution and star forma- 
tion history of the universe (i.e., dn/dt oc 1 It and homo- 
;eneo us thr ough all spa c e, resp ectively), iGal- Yam et al.1 
2005) and iNakar et al.l (|2005l ) have compared whether 
their set of models can produce results statistically con- 
sistent with observations. Among other things they con- 
clude that, within these conventional assumptions, the 
merger model seems inconsistent with the data. 

These previous predictions assume homogeneous star 
forming conditions throughout the universe, with rate 
proportional to the observed time-dep endent star for- 
matio n rate (as given by, for example, Nagamin e et al.l 
(2006) and references therein). In reality, however, 
the universe is markedly heterogeneous as well as time- 
dependent; for example, large elliptical galaxies form 
most of their stars early on. Similarly, predictions for the 
delay time distribution and indeed the total number of 
compact binaries depend strongly on the assumptions en- 
tering into population synthesis simulations. These simu- 
lations evolve a set of representative stellar systems using 
the best parameterized recipies for weakly-constrained 
(i.e., supernovae) or computationally taxing (i.e., stellar 
evolution) features of single and binary stellar evolution. 
By changing the specific assumptions used in these re- 
cipies, physical predictions such as the NS-NS merger 



rate can vary by a few orders of magnitude (see, e.g. 
iKalogera et al.l I2001L and references therein). In partic- 
ular, certain model parameters may allow much better 
agreement with observations. 

In this study we examine predictions based on a large 
database of conventional binary population synthesis 
models: two sets of 500 concrete pairs of simulations 
O, each pair of simulations modeling elliptical and 
spiral galaxies respectively. 1 In order to make predic- 
tions regarding the elliptical-to-spiral rate ratio for bi- 
nary mergers, we adopt a two-component model for the 
star formation history of the universe ($3T]) . Our predic- 
tions include many models which agree with all existing 
(albeit weak) observational constraints we could reliably 
impose. Specifically, many models (roughly half of all 
examined) reproduce the observed short-GRB redshift 
distribution, when we assume either NS-NS or BH-NS 
progenitors. Fewer NS-NS models (roughly a tenth of all 
examined) can reproduce both the short GRB redshift 
distribution and the NS-NS merger rate in spiral-type 
galaxies, as inferred from observ ations of double p ulsars 
seen in the Milky Way (see,e.g. IKim et al.l l2003h . We 
extensively describe the properties of those simulations 
which reproduce observations (@: the redshift distri- 
bution, the fraction of bursts with spiral hosts, and the 
expected detection rate (given a fixed minimum burst 
luminosity). We present our conclusions in section [6] 



2. GAMMA RAY BURSTS: SEARCHES AND 
OBSERVATIONS 

2.1. Emission and detection models 

To compare the predictions of population synthesis cal- 
culations with the observed sample of short GRBs, we 
must estimate the probability of detecting a given burst. 
We therefore introduce (i) a GRB emission model con- 
sisting of an effective luminosity function for the isotropic 
energy emitted, to determine the relative probability 
of various peak fluxes, and a spectral model, for K- 
corrections to observed peak fluxes, and (ii) a detection 
model introducing a fixed peak-flux detection threshold. 
Overall we limit attention to relatively simple models 
for both GRB emission and detection. Specifically, we 
assume telescopes such as BATSE and Swift detect all 
sources in their time-averaged field of view (~ 27T and 
1.4 steradians, respectively; corresponding to a detector- 
orientation correction factor fd given by 1/ fd = 1/2 and 
1.4/4-71") with peak fluxes at the detector Fd greater than 
some fixed t hreshold of Fd = lph sec~ 1 cm~ 2 in 50 to 300 
keV (see.e.g. lHakkila et al"1ll997| ). We note that Swift's 
triggering mechanism is more complex (Gehrels, private 
communication) and appears biased against detections of 
short bursts; for this reason, BATSE results and detec- 
tion strategies will be emphasized heavily in what follows. 

Similarly, though observations of short gamma 
ray bursts reveal a variety of spectra (see, e.g. 
iGhirlanda et al.ll2004L keeping in mind the observed peak 
energy is redshifted), and though this variety can have 

1 Because simulations that produce many BH-NS mergers need 
not produce many NS-NS mergers and vice-versa, we perform two 
independent sets of 500 pairs of simulations, each set designed to 
explore the properties of one particular merger type (i.e. BH-NS or 
NS-NS). The statistical biases motivating this substantial increase 
in computational effort are discussed in the Appendix. 
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Fig. 1. — Characteristic distance to a source y N/4nrFd ver- 
sus its comoving distance. Points: Short bursts with well-defined 
redshifts (SHI; see Table [TJ. Solid line: The critical character- 
istic distance r c (z) = \J N d (z)/4irF d = r{z)yj (1 + z)k(z) versus 
comoving distance r(z), for our simple power-law spectral model 
F v oc v 0,5 . Given our assumptions, systems with fluxes N corre- 
sponding to points above this curve can be seen at the earth with 
a band-limited detector in 50 — 300 keV with peak flux > F d . 



significant implications for the detection of moderate- 
redshift [z > 1) bursts, for the purposes of this paper we 
assume all short gamma ray bursts possess a pure power- 
law spectrum F v oc v~ a with a = —0.5. Th o ugh s everal 
authors such as lAndol (|2004l ) and ISchmidtJ (1200 If ) have 
employed more realistic bounded spectra, similar pure 
power-law spectra have been applied to interpret low- 
redshift o bservations in previ ous theoretical data analy- 
sis effo rts: |Nakar et "aLl (12005T) uses precisely this spectral 
index; iGuetta fc Piranl (|2006l) use a = -0.1. 2 

Because our spectral model is manifestly unphysical 
outside our detection band (50 — 300 keV), we cannot 
relate observed, redshifted fluxes to total luminosity. In- 
stead, we characterize the source's intrinsic photon lumi- 
nosity by the rate N = dN / dt e at which it appears to 
produce B = 50 — 300 keV photons isotropically in its 
rest frame, which we estimate from observed fluxes F in 
this band via a K-correction: 



k(z)i 



:F(47rr 2 )(l + z)k(z) 
! B F v dv/v 



Ib(1+z) F » dv l V 



= (1 + *) 



-0.5 



(1) 

(2) 



To 

-w 



where r(z) is the comoving distance at redshift z. 
give a sense of scale, a luminosity L/(10 erg _1 o 
corresponds to a photon luminosity AT/(4 x 10 53 s _1 ); 
similarly, the characteristic distance to which a photon 
flux can be seen 



is 



1053 s -l)l/2( F(J / lcm 



r c ee ^jN/AnF d ~ 57Mpc(A>/4 x 



-2„-l^-l/2 



Finally, we assume that short GRBs possess an intrin- 
sic power-law peak flux distribution: that the peak fluxes 
seen by detectors placed at a fixed distance but random 
orientation relative to all short GRBs should either (i) be 
precisely zero, with probability 1 — 1/ fb or (ii) collectively 

2 In reality, however, a break in the spectrum is often observed, 
redshifted into the detection band. Under these circumstances, the 
K-correction can play a significant role in detectability. 



be power-law distributed, from some (unknown) mini- 
mum peak flux to infinity, with some probability 1/ fb- 
[This defines fb, the beaming correction factor, in terms 
of the relative probabilities of a visible orientation.] For 
convenience in calculation, we will convert this power- 
law peak-flux distribution into its equivalent power-law 
photon rate TV distribution 



P(> N) = { ^[(N/Nnun) 1 ^ if N > N, 
I fb 



(3) 



if AT < AU, 

where we assume (3 = 2; this particular choice of the 
power-law exponent is a good match to the observed 
BATSE peak-flux distribution (see, e.g. IGuetta &: Piranl 
[200l iNakar et all 120051 : [Aridol l200l ISchmidtJ I2001L and 
references therein). The fraction of short bursts that are 
visible at a redshift z is thus P(z) ee P(> Nj), where Nd 
is shorthand for 4irr 2 (l + z)k(z)Fd- Once again, these as- 
sumptions correspond approximately to those previously 
published in the literature; elementary extensions (for ex- 
ample, a wider range of source luminosity distributions) 
have been successfully applied to match the observed 
BATSE flux distributions and Swift redshift-luminosity 
data [ e.g., in addition to the references mentioned previ- 
ously, IGuetta fc Piranl (|2005t) ] . 

2.2. GRB Observations 

While the above discussion summarizes the most crit- 
ical selection effects - the conditions needed for GRB 
detection - other more subtle selection effects can signif- 
icantly influence data interpretation. Even assigning a 
burst to the "short" class uses a fairly coarse phenomeno- 
logical classification [compare , e.g., the modern spec - 



tral evolution classification of iNorris fc Bonnell 
the machine-learning methods of lHakkila et al 



2006) 



2003), 



and t he original classification paper iKouveliotou et all 
(1993)]; alternate classifications produce slightly but sig- 
nifica ntly different populations (see, e.g. IDonaghv et al.l 
I2006L for a concrete, much broader classification scheme). 
Additionally, short GRB redshift measurements can be 
produced only after a second optical search, with its own 
strong selection bia ses toward low-redshift hosts (see, e.g., 
IBerger et aHl2006D . 

To avoid controversy, we therefore assemble our list 
of short GRBs from four pre viously-published compila- 
tions: (i) IBerger et all (|2006h (Table 1), which provides 
a state-of-the-art Swift-dominated sample with relatively 
homogeneous selection effects; (ii) Donag hv et al.l (|2006f ) 
(Table 8), a broader sample defined using an alte r native 
short-long classification : and finally (iii) IBerger! (|2007l ) 
and (iv) iGehrels et all (|2007f ) which cover the gaps be- 
tween the first two and the present. [We limit attention 
to bursts seen since 2005, so selection effects are fairly 
constant through the observation interval. For similar 
reasons, we do not i nclude the post-fac to IPN galaxy as- 
sociations shown in INakar et all (|2005h (Tab le 1).] This 
compi lation omits GRB 050911 discussed in I Page et al.l 
( 2006) but otherwise includes most proposed short GRB 
candidates. As shown in Table [TJ the sample consists 
of 21 bursts; though most (15) have some redshift in- 
formation, only 11 have relatively well-determined red- 
shifts. However, even among these 12 sources, some dis- 
agreement exists regarding the correct host associations 
and redshifts of G RBs 060505 and 060502B (see,e.g., 
IBerger et aLll2006D . 
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"Gamma-ray burst index 

'Detector in which the GRB was initially detected; S denotes Swift, H denotes HETE-II. 
c Redshift of the host, if well identified. 
^Duration of the burst. 

c Pcak photon flux of the burst (ph/cm 2 /s). 
^Whether the host was optically identified. 
9 Whether the burst produced a visible optical afterglow. 
^Morphology of the host: elliptical (E) or spiral (S). 

•Summary of the previous columns: SI bursts were initially seen by Swift and have a well-defined redshift; S2 bursts were seen by Swift 
and have some uncertain redshift information; S3 bursts include all bursts seen by Swift only. Similarly, SHI includes all bursts seen by 
Swift or HETE-2 with a well-defined red shift. 

^' References: (1 ) IDonaghv eFHTI ||2006T) (2) IGehrels et"al| <2003) f3) ILee et al.l (120051) (41 IBloom eTaT] H2006bl) (5) IBerger et al.l d2005al) 
(6)IBergerl 1120071) fYHBarthelmv et al.l (2005) f8 VFox et ail (120051 ) (91 IVillasenor et al.l (120051 ) (lO )IPedersen et aU 112003) jlD ICovino et al. 
( 200611 (12) IGehrels et alj (I2007D (131 [BeTger et al.l d2005bf) (Ml lProchaska et alJ <200^P (15') ICampana et"aTT (1200611 ( 16) IGrupe et al. 
(2006a) (17) Berger (2006) (18) Nakar (2007) (19) La Parola ct al. (2006) (20) Dietz (2006) (21) Bcrgcr ct al. (2006) (22) Soderbcrg ct al. 
1 2006 jj (2 3) ILevan et al.1 (120061 ) (24) Ide Ugarte Postigo et al.l (120061 ) (251 IRomingl ll2(5oa l (261 IBloom et al.ll2^06al1 (27) IThoene et al. 
1 200711 (28) IQfek et al.l I2007D 



To make sure the many hidden uncertainties and se- 
lection biases are explicitly yet compactly noted in sub- 
sequent calculations, we introduce a simple hierarchical 
classification for bursts seen since 2005: Sn represent the 
bursts detected only with Swift; SHn the bursts seen ei- 
ther by Swift or HETE-II; n = 1 corresponds to bursts 
with well-determined redshifts; n — 2 corrresponds to 
bursts with some strong redshift constraints; and n = 3 
includes all bursts. 

Starting in May 2005, Swift detected 9 short GRBs 
in a calendar year. For the purposes of comparison, 
we will assume the Swift short GRB detection rate to 
be Rd, Swift = 10yr ; compensating for its partial sky 
coverage, this rate corresponds to an all-sky event rate 
at earth of /d,Swifti?D, Swift — 90 yr -1 . However, in or- 
der to better account for the strong selection biases ap- 
parently introduced by the Swift triggering mechanism 
against short GRBs (Gehrels, private communication), 
we assume the rate of GRB events above this threshold 
at earth to be much better represented by the BATSE 
detection rate -R^batse when corrected for detector 
sky coverage, name ly /cz,batse^?d,batse = 170 yr -1 
(|Paciesas et al.lu~999t ) 3 . For similar reasons, in this paper 



we express detection and sensitivity limits in a BATSE 
band (50-300 keV) rather than the Swift BAT band. 

2.3. Cumulative redshift distribution 

As iNakar et al.l (|2005f) demonstrated and as described 
in detail in 21 the cumulative redshift distribution de- 
pends very weakly on most parameters in the short GRB 
emission and detection model (i.e., fa, N, and Fd)- 
When sufficiently many unbiased redshift measurements 
are available to estimate it, the observed redshift distri- 
bution can stringently constrain models which purport 
to reproduce it. At present, however, only 11 reliable 
redshifts are available, leading to the cumulative redshift 
distribution shown in Figure [2] (thick solid line). We 
can increase this sample marginally by including more 
weakly-constrained sources. In Figure [2] (shaded region) 
we show several distributions consistent with SH2, choos- 
ing redshifts uniformly from the intersection of the region 
satisfying any constraints and < z < 5 (an interval 
which encompasses all proposed short GRB redshifts). 
Because this larger sample includes a disproportionate 
number of higher-redshift possibilities, the resulting cu- 



3 Section 2 of Guctta & Piran ( 2005) describes how this rate can 



be extracted from the BATSE catalog paper, taking into account 
time-dependent changes in the instrument's selection effects. 
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Fig. 2. — Cumulative redshift distribution of detected short 
GRBs. The thick solid curve provides the cumulative distribution 
of well-constrained GRBs (i.e., the class SHI). The shaded re- 
gion indicates the range of cumulative distributions produced by 
assigning redshifts to the weakly-constrained (i.e., the class SH2) 
in a manner consistent with the constraints. When only an upper 
or lower limit is available, we pick redshifts using a uniform prior 
for redshifts between and 5. 



mulative redshift distributions still agree at very low red- 
shifts. 

The small sample size seriously limits our ability to ac- 
curately measure the cumulative distribution: given the 
sample size, a Kolmogorov-Smirnov 95% confidence in- 
terval includes any distribution which deviates by less 
than 0.375 from the observed cumulative distribution. 
Rather than account for all possibilities allowed by ob- 
servations, we will accept any model with maximum dis- 
tance less than 0.375 from the cumulative redshift dis- 
tribution for the well-known bursts (i.e., from the solid 
curve in in Figure (H) . 

By performing deep optic al searches to i dentify hosts 
for unconstrained bursts, iBerger et alj (|2006j ) have 
demonstrated that recent afterglow studies are biased to- 
wards low redshift - nearby galaxies are much easier to 
detect optically than high-redshift hosts - and that a sub- 
stantial population of high-redshift short bursts should 
exist. This high-redshift population becomes more ap- 
parent when a few high-redshift afterg lows seen with 
HETE -II before 2005 are included; see Ibonagh v et al.l 
(|2006l ) for details. 

2.4. Comparison with prior work 

Short GRB interpretation: Several previous efforts have 
been made to test quantitative MH-based predictions 
for t h e host, redshif t , lum in osity, and age d i stribu - 
tions [Meszaros et all (120061); iGuetta fc Piranl (120061): 
INakar et alj (12001): IGal-Yam et all (120051): [Bloom et al l 
1999): IBelczvnski et alj (|2006cD : iPerna fc Belczvnskil 
2002)]. However, many authors found puzzling 
discrepanc ies; most not a bly, as has been e mpha - 
sized bv IGal-Yam et al.l (120051) : INakar et alJ ((20051 ) 



and by IGuetta fc Piranl (|2006T ) (bv comparing redshift- 
luminosity distributions to models) and as has seemingly 
been experimentally corroborated with GRB 060502B 
iBloom et al.l (|2006al) . typical observed short GRBs ap- 
pear to occur « (1 — few) x Gyr after their progenitors' 
birth. By contrast, most authors expect population syn- 
thesis predicts a delay time distribution dp/dt cx 1/t 
(e.g., Piran 1992), which was interpreted to imply that 



short delay times dominate, that DCO mergers occur 
very soon after birth, and that mergers observed on our 
light cone predominantly arise from recent star forma- 
tion. Additionally, on the basis of the observed r e dshift- 
lumin osity distri b ution alone, IGuetta fc Piranl (|2006l ) 
and INakar et alj ((2005) conclude short GRB rates to 
be at least comparable to observed present-day NS-NS 
merger rate in the Milky Way. They both also note 
that a large population of low-luminosity bursts (i.e., 
L < 10 49 erg) would remain undetected, a possibility 
which may have some experimental support: post-facto 
correlations between short GRBs and nearby galaxies 
suggests the minimum luminosity of gamm a ray bursts 
(Lmin) could be orders of magnitude lower (|Nakar et al.l 
120051 : iTanvir et al.ll2005l ) . Such a large population would 
lead to a discrepancy between the two types of inferred 
rates. In summary, straightforward studies of the ob- 
served SHB sample suggest (i) delay times and (ii) to 
a lesser extent rate densities are at least marginally 
and possibly significantly incongruent with the observed 
present-day Milky Way sample of double NS binaries, 
and by exte nsion the merger h ypothesis (cf. Sections 
3.2 and 4 of INakar et al.l 12005). A more recent study 
bv lBerger et alj ( 20061 ) suggests that high-redshift hosts 
may be significantly more difficult to identify optically. 
Using the relatively weak constraints they obtain regard- 
ing the redshifts of previously-unidentified bursts, they 
reanalyze the data to find delay time distributions con- 
sistent with dP/dt cx 1/t, as qualitatively expected from 
detailed evolutionary simulations. 

In all cases, however, these comparisons were based 
on elementary, semianalytic population models, with no 
prior on the relative likelihood of different models: a 
model with a Gyr characteristic delay between birth and 
merger was a priori as likely as dP/dt cx 1 jt. For this rea- 
son, our study uses a large array of concrete population 
synthesis simulations, in order to estimate the relative 
likelihood of different delay time distributions. 
Population synthesis: Earlier population synthesis stud- 
ies have explored similar subject matter, even in- 
cluding heterogeneous galaxy populations (see , e.g . 
Belczvnski et al.1 l2006d: Ide Freitas Pacheo et al.l 1200 
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Perna fc B elczvnski 2002; B loom et al.lll999l : lFrver et al 
1999; Bel czvnski et al.ll2002a ). These studies largely ex- 
plored a single preferred model, in order to produce what 
these authors expect as the most likely predictions, such 
as for the offsets expected from merging supernovae- 
kicked binaries and the likely gas content of the circum- 
burst environment. Though preliminary short GRB ob- 
servati ons appeared to con tain an overabundance of short 
GRBs (INakar et al.ll2005D. recent observational analyses 
such as lBerger et al.l (|2006l) suggest high-redshift bursts 
are also present, in qualitative agreem ent with the de- 
tailed population synthesis study by IBelczvnski et al.l 
(2006q). The present study quantitatively reinforces this 
conclusion through carefully reanalyzing the implications 
of short GRB observations, and particularly through 
properly accounting for the small short GRB sample size. 

The extensive parameter study included here, however, 
b ears closest relation to a similar slightly smaller study 
in IBelczvnski et all (|2002af l. based on 30 population syn- 
thesis models. Though intended for all GRBs, the range 
of predictions remain highly pertinent for the short GRB 
population. In most respects this earlier study was much 
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broader than the present work: it examined a much wider 
range of potential central engines (e.g., white dwarf-black 
hole mergers) and extracted a wider range of predictions 
(e.g., offsets from the central host). The present paper, 
however, not only explores a much larger set of popu- 
lation synthesis models (~500) - including an entirely 
new degree of freedom, the relative proportion of short 
GRBs hosted in elliptical and spiral galaxies - but also 
compares predictions specifically against short GRB ob- 
servations. 

3. OTHER RELEVANT OBSERVATIONS 
3.1. Multicomponent star formation history 

The star formation history of the universe has been 
extensively explored through a variety of methods: ex- 
traglactic backg r ound l ight modulo ex tinction (see, e.g., 
iNagamine et al.l 120061 : I Hopkins! 12004 and references 
therein); direct g alaxy counts augm ented by mass 
estimates (see, e.g. iBundv et al.l I2005L and references 
therein); galaxy counts augmented with reconstructed 
star formati on histories from the i r spectral energy distri- 
bution (e.g. iHeayens et al.1 120041: iThompson et all 120061: 
lYan et al.l l2006t lHanish et al.ll2006h: an d via more gen- 
eral composite methods ( Strigari et al.l lifOOS'). Since all 
methods estimate the total mass formed in stars from 
some detectable quantity, the result depends sensitively 
on the assumed low-mass IMF and ofte n on extinction . 
Ho wever, as recently dem onstrated by iHopkinsI (2006) 
and INagamine et al ."(2006). when observations are inter- 
preted in light of a common Chabrier IMF, observations 
largely converge upon a unique star-formation rate per 
unit comoving volume p = dM/dVdt bridging nearby 
and distant universe, as shown in Figure [3] 

Less clearly characterized in the literature are the com- 
ponents of the net star formation history p: the history 
of star formation in relatively well-defined subpopula- 
tions such as elliptical and spiral galaxies. 4 For most 
of time, galaxies have existed in relatively well-defined 
populations, with fairly little morphologi cal evolution 
outside of rare overdense regions (see, e.g. IBundv et al.l 
2005; Ha nish et al.l I2001H and references therein). Dif- 
ferent populations possess dramatically different histo- 
ries: the most massive galaxies form most o f their stars 
very early on (see, e.g. iFeulner et al.l I2005D and hence 
at a characteristically lower metallicity. Further, as 
has been extensively advocated by Kroup a (see, e.g. 
iKroupa & Weidnerl 120031 : iFardal et all 120061 and refer- 
ences therein) the most massive structures could con- 
ceivably form stars through an entirely different collapse 
mechanism ( "starburst-mode" , driven for example by 
galaxy collisions and capture) than the throttled mode 
relevant to disks of spiral galaxies ( "disk- mode" ) , result- 
ing in particular in a different IMF. 

Both components significantly influence the present- 
day merger rate. For example, the initial mass func- 
tion determines how many progenitors of compact bina- 
ries are born from star-forming gas and thus are avail- 
able to evolve into merging BH-NS or NS-NS binaries. 
Specifically, as shown in detail in £14.11 and particularly 

4 Short GRBs have been associated with more refined morpho- 
logical types, such as dwarf irregular galaxies. For the purposes of 
this paper, these galaxies are sufficiently star forming to be "spiral- 
like". 
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Fig. 3. — Star formation history of the universe used in this 
paper. Solid line: Net star formation history implied by Eq. Q. 
Dashed, dotted line: The star formation history due to elliptical 
and spiral galaxies, respectively. 



via Figure 0] elliptical galaxies produce roughly three 
times more high mass binaries per unit mass than their 
spiral counterparts. Add i tional ly, as first recognized by 
Ide Freitas Pacheco et al.l (|2006| ). even though elliptical 
galaxies are quiescent now, the number of compact bina- 
ries formed in ellipticals decays roughly logarithmically 
with time (i.e., dn/dt oc 1/t). Therefore, due to the high 
star formation rate in elliptical-type galaxies ~ 10 Gyr 
ago, the star forming mass density Sp e born in ellipticals 
roughly t e ~ 10 Gyr ago produces mergers at a rate den- 
sity ~ 6p e /t e that is often comparable to or larger than 
the rate density of mergers occurring soon after their 
birth in spiral galaxies ~ dp s /dt. 

Standard two-component model 

As a reference model we use the two -component star 
forma tion history model presented by INagamin e et al] 
(2006). This model consists of an early "elliptical" com- 
ponent and a fairly steady "spiral" component, with star 
formation rates given by 



p = p e + p s 
pc = Ac(t/T C )e- f/TC 



(4) 
(5) 



where cosmological time t is measured starting from 
the beginning of the universe and where the two 
components decay timescales are r e>s — 1.5 and 
4.5 Gyr, respectively (s ee Section 2 and Table 2 
of INagamine et al.l l2006h . These normalization con- 
sta nts A e . s = 0.2 3 , 0.lS M^yr^Mpc" 3 were chosen 
by INagamine et al.l (|2006h so the integrated amount 
of elliptical and spiral star formation agrees with (i) 
the cosmological b a ryon ce nsus [17 * ~ 0.00 2 2; see 
iFukugita fc Peebles! (|2004l ): iRead fc Trenthaml (|2005l ) 
and references therein]; (ii) the expected degree of gas 
recycling from one generation of stars to the next; and 
(iii) the relative masses in different morphological com- 
ponents (60% : 40%). Explicitly, these two constants 
are chosen so J p e /p c = £7»/0.55 x 0.6 and J p e /p c — 
f2*/0.55 x 0.4, respectively. 

Each component forms stars in its own distinctive con- 
ditions, set by comparison with observations of the Milky 
Way and elliptical galaxies. We assume mass converted 
into stars in the fairly steady "spiral" component is done 
so using solar metallicity and with a fixed high-mass IMF 
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power la w \p = —2.7 in the brok en-power-law Kroupa 
IMF; see lKroupa fc Weidneri (|2003l )]. On the other hand, 
we assume stars born in the "elliptical" component are 
drawn from a broken power law IMF with high-mass in- 
dex within p E [—2.27, —2.06] and metallicity Z within 
0.56 < Z/Zq < 2. These elliptical birth conditions 
agree with obse rvations of bot h old ellipticals in the lo- 
cal universe (see lLi et al.ll2006l and refer ences therein) as 
well as of young st arburst clusters (see iFall et alj [20051 : 
I Zhang fc Fall|[T999l and references therein). 

3.2. Binary pulsar merger rates in the MilkyWay 

If binary neutron stars are the source of short GRBs, 
then the number of short GRBs seen in spirals should 
be intimately connected to the number of binary pul- 
sars in the Milky Way that are close enough to merge 
through the emission of gravitational radiation. Four 
unambiguously merging double pulsars have been found 
within the Milky Way usin g pulsar surveys with well- 
understood selection effects. iKim et aTl (|2003| ) developed 
a statistical method to estimate the likelihood of double 
neutron star formation rate estimates, designed to ac- 
count for the small num ber of known sy s tems and their 
associated uncertainties. iKaTo" sera et~aTl (|200l summa- 
rize the latest results of this analysis: taking into account 
the latest pulsar data, a standard pulsar beaming correc- 
tion factor fi, = 6 for the unknown beaming geometry of 
PSR J0737-3037, and a likely distribution of pulsars in 
the galaxy (their model 6), they constrain the rate to 
be between r MW = 16.9Myr _1 and 292.1Myr~ 1 . (95% 
confidence interval) 5 . 

Assuming all spiral galaxies to form stars similarly to 
our Milky Way, then the merger rate density in spirals 
at present 7?.[ s ] (T) must agree with the product of the 
formation rate per galaxy tmw and the density of spiral 
galaxies n s . Based on the ratio of the blue light den- 
sity of the universe to the blue light attributable to the 
Milky Way, the density of Milky Way- equivalent galaxies 
lies betw e en 0.7 5 x 10~ 2 Mpc~ 3 and 2 x 10~ 2 Mpc~ 3 (see 



iPhinneyl (119911). iKalogera et all ([20011) . iNutzman et all 



( 2004h . iKopparapu et al.l ( 2007T ) and references therein). 
We therefore expect the merger rate density due to 
spirals at present to lie between 0.15Myr _1 Mpc~ 3 and 
5.8Myr _1 Mpc" 3 (with better than 95% confidence). 

4. PREDICTIONS FOR SHORT GRBS 

4.1. Population synthesis simulations 

We study the formation of compact objects with the 
StarTrack popul a tion sy nthesis code, first developed by 
IBelczvnski et al.l (|2002bl ) and recently significantly ex - 
tended as described in detail in lBelczynski et aTl (|2006af ); 

5 The range of binary neutron star merger rates that we ex- 
pect to contains the true present-day rate has continued to evolve 
as our knowledge about existing binary pulsars and the distri- 
bution of pulsars in the galaxy has changed. The range quoted 
here reflects the recent calculations on binary pulsar merger rates, 
and corresponds to t h e merg er rate confidence interval quoted in 
lO'Shau ghncssv et al. ( 20071) ) (albeit with a different convention 
for assigning upper and lower confidence interval boundaries). In 
particular, this estimate does not incorportate conjectures regard- 
ing a possibly sh orter lifetime of PSR J0737-3037, as described in 
IKim et alj :f?IK)(i :. The properties of this pulsar effectively deter- 
mine the present-day merger rate, and small changes in our under- 
standing of those properties can significantly change the confidence 
interval presented. 



see §2 of IBelczvnski et aTl ([2006b) for a succinct descrip- 
tion of the changes between versions. 

Since our understanding of the evolution of single and 
binary stars is incomplete, this code parameterizes sev- 
eral critical physical processes with a great many param- 
eters (~ 30), many of which influence compact-object 
formation dramatically; this is most typical with all 
current binary population synthesis codes used by var- 
ious groups. For the StarTrack population synthesis 
code, in addition to the IMF and metallicity (which vary 
depending on whether a binary is born in an ellipti- 
cal or spiral galaxy), seven parameters strongly influ- 
ence compact object merger rates: the supernova kick 
distribution (modeled as the superposition of two in- 
dependent Maxwellians, using three parameters: one 
parameter for the probability of drawing from each 
Maxwellian, and one to characterize the dispersion of 
each Maxwellian), the solar wind strength, the common- 
envelope energy transfer efficiency, the fraction of an- 
gular momentum lost to infinity in phases of non- 
conservative mass transfer, and the relative distribution 
of masses in the binary. Other parameters, such as 
the fraction of stellar systems which are binary (here, 
we assume all are, i.e., the binary fraction is equal 
to 1) and the distribution of initial binary parameters , 
are comparatively we ll -deter mined (see e.g.|Abt (1983), 
iDuquennov fc Mavorl (|1991h and references therein). 6 
Even for the less-well-constrained parameters, some in- 
ferences have been drawn from previous studies, either 
more or less directly (e.g., via observations of pulsar 
proper motions, which presumably relate closel y to su - 
pernovae kick stren g th; se e, e.g., iHobbs et all (12 005). 
lArzoumanian et aTl (|2002l ). iFaucher-Giguere fc Kaspil 
(|2006h and references therein) or via comparison of 
some subset of binary popul ation synthesis result s 
with observations ( e. g., §8 of IBelczvnski et al.l (|2006al) , 
van der Sluvs et al.l (120061) . iNe lcmans fc ToutJ (l2005li . 



Willems fc Kolbl (l2002D . lPodsiadlowski et all (|2002F and 
references therein). Based on these and other com- 
parisons, while other parameters entering into popula- 
tion synthesis models can influence their results equally 
strongly, these particular seven parameters are the 
least constrained observationally. For this reason, de- 
spite observational suggestions that point towards pre- 
ferred values for these seven parameters - and despite 
the good agreement with short GRB and other ob- 
servations obtain e d whe n using these preferred values 
(|Belczvnski et al.l (|2006d )) - in this paper we will con- 
servatively examine the implications of a plausible range 
of each of these parameters. Mor e specifically, despite 
the Milky W ay-specific studies of lO'Shaughncss v et alJ 
(|2005l 120075 ) (which apply only to spirals, not the ellip- 
tical galaxies included in this paper), in this study we will 
continue to assume all seven parameters are unknown, 
drawn from the plaus i ble par ameter ranges described in 
lO'Shaughnessv et all (|2007bD . 

As noted previously in § 13. 1[ we perform simulations of 
two different classes of star-forming conditions: "spiral" 
conditions, with Z — Zq and a high-mass IMF slope 

6 Particularly for the application at hand - the gravitational- 
wave-dominated delay between binary birth and merger - the de- 
tails of the semimajor axis distribution matter lit t le. For a similar 
but more extreme case, sec O'Shaughncssy ct al. (2007c). 
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of p = —2.7, and "elliptical" conditions, with a much 
flatter IMF slope and a range of allowed metallicities 
0.56 < Z/Zq < 2. 

Archive selection: Our collection of population synthesis 
simulations consists of roughly 3000 and 800 simulations 
under spiral and elliptical conditions, respectively. Our 
archives are highly heterogeneous, with binary sample 
sizes N that spread over a large range. 7 A significant 
fraction of the smaller simulations were run with pa- 
rameters corresponding to low merger rates, and have 
either no BH-NS or no NS-NS merger events. There- 
fore, though the set of all population synthesis simula- 
tions is unbiased, with each member having randomly 
distributed model parameters, the set of all simulations 
with one or more events is slightly biased towards simu- 
lations with higher-than-average merger rates. Further, 
the set of simulations with many events, whose proper- 
ties (such as the merger rate) can be very accurately es- 
timated, can be very strongly biased towards those mod- 
els with high merger rates. Fortunately, as discussed at 
much greater length in the Appendix, the set of simula- 
tions with nN > 5 x 10 6 and n > 20 has small selection 
bias and enough simulations (976 and 737 simulations 
NS-NS and BH-NS binaries under spiral-like conditions, 
as well as 734 and 650 simulations under elliptical condi- 
tions, respectively) to explore the full range of population 
synthesis results, while simultaneously insuring each sim- 
ulation has enough events to allow us to reliably extract 
its results. 

4.2. Results of simulations 

From each population synthesis calculation (a) per- 
formed under elliptical or spiral conditions (C = e, s) 
and for each final result (K), we can estimate: (i) the 
number of final K events per unit mass of binary merg- 
ers progenitors, i.e., the mass efficiency (\c,a,K)\ and 
(ii) the probability P c , a ,K(< t) that given a progenitor 
of K the event K (e.g., a BH-BH merger) occurs by time 
t since the formation of K. Roughly speaking, for each 
simulation we take the observed sample of n binary pro- 
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genitors of K, with Mi., 
estimate 

A = 



and delay times tx...n> and 



n fcut 



' N (M) 

p m (<t)=X;e(t-t J ) 



(6) 
(7) 



where Q(x) is a unit step function; N is the total num- 
ber of binaries simulated, from which the n progenitors 
of K were drawn; (M) is the average mass of all possible 
binary progenitors; and f cut is a correction factor ac- 
counting for the great many very low mass binaries (i.e., 
with primary mass mi < m c = 4M Q ) not included in our 
simulations at all. Expressions for both (M) and f cu t in 
terms of population sy nthesis model parame t ers are pro- 
vided in Eqs. (1-2) of lO'Shaughnessv et all (|2007af) . In 
practice, P m (t) and dP m /dt are estimated with smooth- 
ing kernels, as discussed in Appendix IB1 Given the char- 
acteristic sample sizes involved (e.g., n > 200 for NS-NS), 

7 In practice, the sample size is often chosen to insure a fixed 
number of some type of event. As a result, usually the sample size 
N and the number of any type of event n are correlated. 
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Fig. 4. — Smoothed histograms of the mass efficiency A [Eq. (0] 
of simulations used in our calculations, shown for spiral (solid) 
and elliptical (dotted) birth conditions. As expected given the 
differences in the IMF, elliptical galaxies produce BH-NS binaries 
roughly three times more efficiently than spirals. However, appar- 
ently because our population synthesis sample involves high ly cor - 
rclatcd choices for n and N (see the Appendix and Figure I A 131 ) , 
our distribution of NS-NS mass efficiencies remains biased, pro- 
ducing identical distributions for both elliptical and spiral birth 
conditions. 



we expect P m to have absolute error less than 0.05 at each 
point (95% confidence) and dP m /dt to have rms relative 
error less than 20% (95% confidence). Since these errors 
are very small in comparison to uncertainties in other 
quantities in our final calculations (e.g., the star forma- 
tion rate of the universe), we henceforth ignore errors in 
P and dP/dt. 

Figures |4l and [6] show explicit results drawn from 
these calculations. From these figures, we draw the fol- 
lowing conclusions: 

Uncertainties in binary evolution significantly affect re- 
sults: As clearly seen by the range of possiblities allowed 
in Figures [3] and [5l our imperfect understanding of binary 
evolution implies we must permit and consider models 
with a wide range of mass efficiencies A and delay time 
distributions P m (< t). 

The merger time distribution is often well-approximated 
with a one parameter family of distributions, dP m /dt oc 
1/t: As suggested by the abundance of near-linear dis- 
tributions in Figure O the delay time distribution P m is 
almost always linear log£. Further, from the relative in- 
frequency of curve crossings, the slope of P m versus logt 
seems nearly constant. As shown in the bottom panels of 
Figure O this constant slope shows up as a strong corre- 
lation between the times £(5%) and £(50%) at which P m 
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reaches 0.05 and 0.5 when logt(5%) /Myr > 1.5: 

W^5n%l ~ / l0 S*( 5% ) + 2 - 5 if logt(5%) > 1.5 
logt^u/oj ~ | ioi gt(5%) - 11 if logt(5%) < 1.5 

(8) 

The merger time distribution is at most weakly correlated 
with the mass efficiency: Finally, as seen in the top pan- 
els of Figure a wide range of efficiencies are consistent 
with each delay time distribution. The maximum and 
minmum mass efficiency permitted increase marginally 
with longer median delay times £(50%) - roughly an or- 
der of magnitude over five orders of magnitude of t(50%). 
But to a good approximation, the mass efficiencies and 
delay times seem to be uncorrelated. 

4.3. Converting simulations to predictions 

Each combination of heterogeneous population synthe- 
sis model, star formation history, and source model leads 
to a specific physical signature, embedded in observables 
such as the relative frequencies with which short GRBs 
appear in elliptical and spirals, the average age of progen- 
itors in each component, and the observed distribution 
of sources' redshifts and luminosities. All of these quan- 
tities, in turn, follow from the two functions 1Zc(t), the 
merger rate per unit comoving volume in C = elliptical 
or spiral galaxies. 

The rate of merger events per unit comoving volume is 
uniquely determined from (i) the SFR of the components 
of the universe dp{c}/dt; (ii) the mass efficiency Ajq at 
which K mergers occur in each component C; and (iii) 
the probability distribution dP m {c}/dt for merger events 
to occur after a delay time t after star formation: 



[C] 



(9) 



dtb^{C,K}- 



dP, 



"' !,7} (*-* 6 )^(4)(io) 



dt y ' dt 
Though usually H]nAt) is experimentally inaccessible, 
because our source and detection models treat elliptical 
and spiral hosts identically, 8 the ratio uniquely deter- 
mines the fraction f s of spiral hosts at a given redshift: 

f s (z)=TZ [s] /(TZ le] +n ls] ) (11) 

Additionally, as described in £13.21 observations of NS- 
NS in the Milky Way constrain the present-day merger 
rate of short GRB progenitors (7?.[ s ] (T un i verse )), if those 
progenitors are double neutron star binaries. 

Unfortunately, the relatively well-understood physical 
merger rate distributions R{c\ are disguised by strong 
observational selection effects described in § [21 notably 
in the short GRB luminosity function. Based on our 
canonical source model, we predict the detection rate 
Rd of short GRBs from to be given by 

(12) 



Rp = / J Rd[c] 



R D [C]=fd / n [c] P [c] (z)4nr 2 cdt 



(13) 



N„ 



fdfbFd 



cdt ■ 



(l + z)k(z) 



8 In practice, gas-poor elliptical hosts should produce weaker 
afterglows. Since afterglows are essential for host identification and 
redshifts, elliptical hosts may be under-represented in the observed 
sample. 



where the latter approximation holds for reasonable 
N m in/ fb < 10 57 s _1 (i.e., below values corresponding to 
observed short bursts). While the detection rate depends 
sensitively on the source and detector models, within the 
context of our source model the differential redshift dis- 
tribution p{z) 



P{) dz^(l + z)k(z)f b 



(14) 



and the cumulative redshift distribution P(< z) = 
f n p{z)dz do not depend on the source or detector model 
(|Nakar et alJl2005h . 

Detected luminosity distribution: In order to be self- 
consistent, the predicted luminosity distribution should 
agree with the observed peak-flux distribution seen by 
BATSE. However, so long as N m in is small for all popula- 
tions, our hypothesized power-law form N~ 2 for the GRB 
luminosity function immediately implies the detected cu- 
mulative peak flux distribution P(> F) is precisely con- 
sistent: P(> F) = (Fd/F), independent of source pop- 
ulatio n; see for example the discussion in iNakar et al.l 

(2005) near their Eq. (5). While more general source 
population models produce observed luminosity func- 
tions that contain some information a bout the redshift 
distrib ution of bursts - for example, iGuetta &: Piranl 

(2006) and references therein employ broken power-law 
luminosity distributions; alternatively, models could in- 
troduce correlations between brightness and spectrum - 
so long as most sources remain unresolved (i.e., small 
N m in), the observed peak flux distribution largely reflects 
the intrinsic brightness distribution of sources. Since 
INakar et~a l. (2005) demonstrated this particular bright- 
ness distribution corresponds to the observed BATSE 
flux distribution, we learn nothing new by comparing 
the observed peak flux distribution with observations and 
henceforth omit it. 

4.4. Predictions for short GRBs 

Given two sets of population synthesis model param- 
eters - each of which describe star formation in ellipti- 
cal and spiral galaxies, respectively - the tools described 
above provide a way to extract merger and GRB de- 
tection rates, assuming all BH-NS or all NS-NS mergers 
produce (possibly undetected) short GRB events. Rather 
than explore the (fifteen-dimensional: seven parameters 
for spirals and eight, including mctallicity, for ellipticals) 
model space exhaustively, we explore a limited array 
of 500 "model universes" by (i) randomly 10 selecting 

9 We draw our two-simulation "model universe" from two collec- 
tions of order 800 simulations that satisfy the constraints described 
in the Appendix, one for ellipticals and one for spirals. Computa- 
tional inefficiencies in our postprocessing pipeline prevent us from 
performing a thorough survey of all ~ 10 5 possible combinations 
of existing spiral and elliptical simulations. 

10 At present, our population synthesis archives for elliptical and 
spiral populations were largely distributed independently; we can- 
not choose pairs of models with similar or identical parameters for, 
say, supernovae kick strength distributions. The results of binary 
evolution from elliptical and spiral star forming conditions, if any- 
thing, could be substantially more correlated than we can presently 
model. We note however that there is no a priori expectation nor 
evidence that evolutionary parameters should indeed be similar in 
elliptical and spiral galaxies. 
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Fig. 5. — Cumulative probabilities P m (< t) that a NS-NS binary (left panels) or BH-NS binary (right panels) will merge in time less than 
t, for twenty randomly-chosen population synthesis models, given spiral (top) and elliptical (bottom) star forming conditions. A vertical 
dashed line indicates the age of the universe. For the sample sizes involved, these distributions are on average accurate to within 0.05 
everywhere (with 90% probability); see Figure [B14I 
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Fig. 6. — Scatter plots to indicate correlations between results of various simulations. Top panels: Scatter plot of mass efficiency log A 
and average delay time logt(50%). Bottom panels: Scatter plot of logt(5%) versus logi(50%); also shown is the analytic estimate of Eq. 
J5). Left panels indicate spiral conditions; right panels indicate elliptical conditions. Despite the differences between these two types of 
simulations (i.e., metallicity and initial mass function), the range of delay time distributions and mass efficiencies are largely similar (i.e., 
the corresponding left and right panels resemble one another). 



two population synthesis simulations, one each associ- 
ated with elliptical (e) and spiral (s) conditions; (ii) esti- 
mating for each simulation the mass efficiency (X e , s ) and 
delay time distributions (P e ,s)\ (hi) constructing the net 
merger rate distribution 1Z using Eqs. (|9I10|) : and (iv) 
integrating out the associated expected redshift distribu- 
tion p(z) [Eq. CL1]. 

The results of these calculations are presented in Fig- 
ures [3 [3 [U [lOl and QT] [These figures also compare 
our calculations' range of results to observations of short 
GRBs (summa rized in Table [H and merging Milky Way 
binary pulsars (|Kim et al.ll2003l ); these comparisons will 
be discussed extensively in the next section.] More 



specifically, these five figures illustrate the distribution of 
the following four quantities that we extract from each 
"model universe" : 

Binary merger rates in present-day spiral galaxies: To 
enable a clear comparison between our multi-component 
calculations, which include both spiral and elliptical 
galaxies, and other merger rate estimates that incorpo- 
rate only present-day star formation in Milky Way-like 
galaxies, the two solid curves in Figure [7] show the distri- 
butions of present-day NS-NS and BH-NS merger rates 
in spiral galaxies seen in the respective sets of 500 simu- 
lations. 

In principle, the BH-NS and NS-NS merger rates 
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Fig. 7. — The distribution of merger rate densities in spiral- type 
galaxies H a for BH-NS mergers (top) and NS-NS mergers (bot- 
tom); the solid curve includes all simulations, the dashed curve 
only those simulations reproducing the observed redshift distribu- 
tion; and the dotted curve (bottom panel only) only those simula- 
tions reproducing the NS-NS merger rate in spiral galaxies derived 
from an anal ysis of binary pulsars (also shown on the bottom panel, 
in gray; see ^13.21 for details); and the dot-dashed curve (top panel 
only) includes only those simulations which, under the most opti- 
mistic assumptions, predict short GRBs should occur at least as 
frequently as has been seen. The bottom pan el in particular should 
be com pared with Figure 3 (top panel) of O'Shauehncssv ct al. 
(20079). 



should be weakly correlated, as the processes (e.g., com- 
mon envelope) which drive double neutron star to merge 
also act on low-mass BH-NS binaries, albeit not always 
similarly; as a trivial example, mass transfer processes 
that force binaries together more efficiently may deplete 
the population of NS-NS binaries in earlier evolution- 
ary phases while simultaneously bringing distant BH-NS 
binaries close enough to merge through gravitational ra- 
diation. Thus, a simulation which contains enough BH- 
NS binaries for us to estimate its delay time distribution 
dP/dt need not have produced similarly many NS-NS bi- 
naries. For this reason, we constructed two independent 
sets of 500 "model universes" , one each for BH-NS or 
NS-NS models. However, as a corollary, the randomly- 
chosen simulations used to construct any given BH-NS 
"model universe" need not have enough merging NS-NS 
to enable us to calculate the present-day merger rate, 
and vice-versa. In particular, we never calculate the 
double neutron star merger rates in the BH-NS model 
universe. Thus, though the BH-NS and NS-NS merger 
rates should exhibit some correlation, we do not explore 
it here. In particular, in the next section where we com- 
pare predictions against observations, we do not require 
the BH-NS "model universes" reproduce the present-day 
NS-NS merger rate. 



Short GRB detection rate: As described in detail in § [2J 
the fraction of short GRBs on our past light cone which 
are not seen depends strongly but simply on unknown 
factors such as the fraction of bursts pointing towards 
us, which we characterize by l//b where fb is the beam- 
ing factor, and the fraction of bursts luminous enough 
to be seen at a given distance, which we characterize by 
P(> N d ) where N d = 4nr 2 k(z)(l + z)F d is the mini- 
mum photon luminosity visible at redshift z. The short 
GRB detection rate also depends on the detector, includ- 
ing the fraction of sky it covers (1//<j) and of course 
the minimum flux F d to which each detector is sensi- 
tive. To remove these significant ambiguities, in Figure 
[5] we use solid curves to plot the distribution of detection 
rates found for each of our 500 "model universes" (top 
panel and bottom panel correspond to the BH-NS and 
NS-NS model universes, respectively), assuming (i) that 
no burst is less luminous than the least luminous burst 
seen, namely, GRB 060505, with an apparent (band- 
limited) isotropic luminosity iV m i nseen ~ 3 x 10 55 s _1 , or 
L~ 7 x 10 48 ergs -1 (see Table [T]); (ii) that beaming has 
been "corrected" , effectively corresponding to assuming 
isotropic detector sensitivity and source emission; and 
(iii) that the detector has a peak flux detection thresh- 
old of F d = lcm _2 s , corresponding roughly to the 
true BATSE and Swift peak flux thresholds presented in 
§ [21 These choices are conservative; therefore, our rate 
estimate for each model universe is an upper bound. 
Short GRB redshift distribution: The short GRB detec- 
tion rate depends on the ability of gamma-ray detec- 
tors to distinguish burst events from background noise 
and the GRB luminosity function. Its selection effects 
therefore depend only on the properties of gamma-ray 
telescopes. On the other hand, the short GRB redshift 
distribution implicitly contains several other selection ef- 
fects which have not been factored into our analysis (e.g., 
regarding the availability and reliability of an association 
of each burst to a host and the ability of optical telescope 
to extract a spectrum and redshift of that host). Lacking 
the ability to characterize these selection effects, we per- 
form the most straightforward prediction for the distribu- 
tion of short GRB redshifts: using Eq. (|14[) . we assume 
all detected short GRB events have accurate and un- 
ambiguous optical follow-up and redshifts. Based upon 
that assumption, we generate two sets of 500 candidate 
redshift distributions for short GRBs, assuming either 
(i) the set of merging BH-NS binaries or (ii) the set of 
merging NS-NS binaries is identical to (iii) the set of all 
short GRBs. The results are shown in the top left panels 
of Figures E] (for NS-NS binaries) and [TO] (for BH-NS bi- 
naries). Rather than provide all 500 cumulative redshift 
distributions, we have sorted these distributions by their 
value at z — 0.1 and plotted indicative curves at chosen 
percentiles: for example, the dotted curves show the 50th 
and 450th cumulative redshift distributions we obtained, 
corresponding to the 10th and 90th percentile. Because 
of the limited number of "model universes" presently 
available to us, we do not presently attempt to resolve 
the tails of the distribution; however the range of po- 
tential cumulative redshift distributions includes a low a 
priori probability tail with members which are slightly 
more biased towards low redshift than the cumulative 
distributions those shown here. 
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Fig. 8. — Distribution of predicted all-sky, beaming- corrected 
short GRB detection rates log Rofbfdi if bursts arise from BH-NS 
(top) and NS-NS (bottom) mergers and all bursts produce a higher 
peak flux of 50-300 keV photons than the least-luminous burst in 
our sample (with TV ~ 3 X 10 55 s" 1 ), compared to the observed 
BATSE (dashed vertical line) and Swift (dotted) all-sky detection 
rate. In other words, a plot of the most optimistic predictions 
of each two-component "model universe" for the burst detection 
rate: assuming a detector has all sky visibility, that all burst emit 
isotropically with no beaming, and that no bursts arc faint enough 
to be missed. Almost all models produce at least as many bursts 
as are observed; all models can therefore reproduce observations by 
choosing some low minimum peak photon luminosity iV m ; n , typi- 
cally 10 -10 3 times lower than the one of the faintest burst seen. As 
in Figure [7] the solid curves include all 500 "model universes"; the 
dashed curves include only those "model universes" which repro- 
duce the short GRB redshift distribution; the dotted curve (bottom 
panel only) includes only models which are reasonably consistent 
with the present-day double neutron star merger rate in the Milky 
Way; and the dot-dashed curve (top panel only) includes only those 
simulations which, under the most optimistic assumptions, predict 
short GRBs should occur at least as frequently as has been seen. 



Fraction of mergers in spiral galaxies: Finally, in Figure 
1111 (solid curves) we show the a priori probability that 
a fraction f s of binary mergers occur in spiral galaxies, 
using the scaled variable X = arctanh(2/ s — 1). This 
representation allows us to better illustrate relative prob- 
abilities of spiral fractions f s that are very near 1 or 0; 
to give a sense of scale, X = 1(2) corresponds to 88% 
(98%) of all mergers occurring in spiral galaxies. 

As seen in Figure 111! a priori we cannot say whether 
elliptical or spiral galaxies should host most presently- 
occurring binary mergers. Furthermore, a significant 
fraction of models have \X\ > 1 (i.e., in those models, 
almost all mergers occur in spirals or ellipticals). This 



extreme range for f s can be traced back directly to the 
large range of mass efficiencies shown in Figure |4] For 
a randomly chosen pair of population synthesis model 
for elliptical and spiral galaxies, one will often be signif- 
icantly larger than the other, leading to a spiral fraction 
f s near its limits (i.e., f s ~ 0, 1). 

5. TESTING PREDICTIONS AGAINST OBSERVATIONS 

To summarize, in this paper we are trying to examine 
whether either of two very simple hypotheses are consis- 
tent with what we have heretofore seen of short GRBs 
and binary mergers are consistent with our understand- 
ing of the star formation history of the universe and of 
binary stellar evolution. The two simplest hypothesis 
are, on the one hand, that every short GRB is a binary 
BH-NS merger and that every BH-NS merger produces 
a short GRB; and on the other hand the correspond- 
ing statement for NS-NS binaries. 11 Further, assuming 
that in both cases some fraction of all possible models 
are consistent with observations, we want to understand 
the properties of those consistent models and the physi- 
cal processes that require those properties to be as they 
are. Of course, these comparisons a re affected by obser- 
vational selection effects. Notably, iBerger et al.l (|2006t ) 
has suggested that the lack of deep optical follow-up on 
faint short GRBs has biased the redshift distribution 
to low redshift. Lacking any ability to control an un- 
known bias, we proceed with the simplest possible test: 
we compare existing results to currently known obser- 
vations. Specifically, we require our predictions for the 
binary neutron star merger rate (Figure [71 available only 
for binary neutron star short GRB "model universes"), 
for the short GRB detection rate (Figure [8|), and for the 
short GRB redshift distribution (Figures [9] and [TOf be 
statistically consistent with the corresponding observa- 
tions, all of which are shown in their respective plots. 
We then compare the distribution of constrained quan- 
tities (in the aforementioned Figures and in Figures 1111 
and[T2"|) with their initial distributions, to understand the 
physics implied by observations. 12 

5.1. Comparisons I: NS-NS as burst source 

No definitive evidence exists to determine the fraction 
of short bursts that are less luminous than the faintest 
currently known; to discover the degree to which, if any, 
short GRBs are beamed 13 ; and to demonstrate that all 
mergers must produce bursts. Therefore, we must take 

11 We start with the simplest and most tractable pair of hy- 
potheses. A more realistic model would allow for a fraction of 
both BH-NS and NS-NS mergers to produce short GRBs, along 
with some proportion of young neutron stars (in bursts from soft 
gamma repeaters, commonly abbreviated SGRs). However, with 
so many degrees of freedom, such a model would be very difficult to 
constrain without additional observational inputs (e.g., direct con- 
firmation of the fraction of mergers which are SGRs in the nearby 
universe; a reliable measure of the fraction of mergers that occur in 
elliptical and spiral galaxies; gravitational wave detection of merger 
events; etc.). 

12 Because of the limited number of constrained models and the 
large number of underlying parameters (fifteen), we do not attempt 
to characterize the underlying physical parameters of constraint- 
satisfying models. Instead, we focus on describing the essential 
features of models, such as a long characteristic delay time or high 
mass efficiency, that allow a model to satisfy observational con- 
straints. 

13 Evidence has been presented to suggest that short GRB emis- 
sion has a break in one or more frequencies, a result that has been 
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Fig. 9. — Demonstration that popul atio n synthesis models can reproduce short GRB redshift distributions, assuming NS-NS mergers 
produce all short GRBs. As in Figure 1101 the top right and left panels compare the range of short GRB redshift distributions expected 
from our two-component model : l%/99% (dotted), 10%/90% (dashed), and 25%/75% (solid) redshift distributions are overlaid on the 
observed cumulative redshift distribution, both for all simulated models (top left) and those models which remain everywhere close to the 
observed redshift distribution (bottom left). Top right: As above, but including only those NS-NS models which produce a NS-NS merger 
rate in agreement with observations of the Milky Way merger rate (Figure . Bottom right: As above, but including only those NS-NS 
models which reproduce the number of binary pulsars seen in the Milky Way and the short GRB redshift distribution. 



the bottom panel of Figure [5] at face value: if no bursts 
are less luminous than those seen, every short GRB 
model based on NS-NS produces many more short GRBs 
than are seen, and therefore every double neutron star 
"model universe" can reproduce the present-day short 
GRB detection rate by a proper choice of, for example, 
the minimum luminosity of short GRBs. Therefore, only 
two of the three available observations - the short GRB 
redshift distribution and the set of merging Milky Way 
binary pulsars - can reject some of the 500 double neu- 
tron star "model universes." 

The statistics of this comparison are summarized in 
Figure [D Relatively few of our double neutron star 
models are consistent with observations of merging dou- 
ble pulsars in the Milky Way 14 : 78 out of 500 models, 
or ~ 16%; see the top right panel of Figure [9] as well 
as Figure [7] Among the few models which agree with 
Milky Way observations, however, most have a cumula- 
tive redshift distribution within 0.375 of our family of 
reference curves and therefore cannot be rejected by a 
Komolgorov-Smirnov test (95% confidence); see the bot- 

interpreted as a "jet break" produced by a beamed jet (see, e.g., 
ISoderbere etld1l2006bl ; IGrupe et al~ll2006bl ; INakarl |2007| , and ref- 
erences therein). At present we choose to remain conservative 
regarding beaming until several multi-band observations confirm 
these breaks exist. 

14 Fewer still would be consistent should we require agreement 
with more binary pu l sar ra te c onstraints, as has been show n in 
IQ'Shaughnessv et all d2005l> and lO'Shaughnessv et alj H2007bl 1. 



torn left panel of Figure [H Therefore 35 out of 500, 
roughly 7% of all models, appear fully consistent with 
all observations considered here and with the assumption 
that short GRBs are produced exclusively by all NS-NS 
mergers. 

Because observations of the Milky Way suggest a NS- 
NS merger rate towards the high end of what our sim- 
ulations produce (see the bottom panel of Fig. [7|), and 
because spiral star formation extends through the recent 
universe - that is, precisely into the time intervals during 
which a significant fraction of short GRBs have been ob- 
served (see Figure [2]) - to an excellent approximation 
we can conclude that, if short GRBs are due to NS- 
NS mergers, the best-fitting models of all observations 
produce mergers and short GRBs preferentially in spiral 
galaxies, with rapid mergers following quickly upon re- 
cent star formation. More specifically, we can gcncrically 
draw the following explicit conclusions: 
High spiral mass efficiency needed: As indicated in Fig- 
ure [71 only a few populaton synthesis simulations can 
produce as many merging NS-NS binaries in spiral galax- 
ies as does a natural extrapolation of the known Milky 
Way NS-NS merger rate. The few "model universes" 
which reproduce these high merger rates necessarily have 
non-typical high mass efficiencies A s for forming double 
neutron stars in spiral galaxies. 

High spiral fraction preferred: As a consequence of such 
high mass efficiencies A s - that is, because of the high 
rate at which the spiral galaxies in these "universes" 
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Fig. 10. — Demonstration that population synthesis models can reproduce short GRB redshift distributions, assuming BH-NS mergers 
produce all short GRBs. The jagged curve and shaded regions provide the cumulative redshift distribution for observed short GRBs (Fig. 
|2f. Top left panel: The smooth curves illustrate the range of short GRB redshift distributions; out of the 500 simulated, the two solid 
curves correspond to the curves with the 25th and 75th percentile values of P(0.1); the dashed curves to the 10th and 90th percentiles; and 
the dotted curves to the first and 99th percentiles. Bottom left panel: As above, but including only those models which produce redshift 
distributions which differ from observations by less than 0.38 (i.e., a 5% Kolmogorov-Smirnov false alarm probability). Top right panel: As 
previously, but including only those model which, when given the most favorable assumptions, still predict too few short GRB detections. 
Bottom right panel: As previously, but requiring both the predicted short GRB rate and redshift distribution be simultaneously consistent 
with GRB observations. 



produce NS-NS systems - most mergers occur in spirals; 
equivalently, the spiral fraction f s of constraint-satisfying 
models is often strongly biased towards spiral galaxies 
(f a > 1/2). The bottom panel of Figure [TT1 explicitly 
shows that, though unconstrained "universes" can pro- 
duce mergers preferentially in either elliptical or spiral 
galaxies (solid curve), in order to reproduce observations 
of the Milky Way's merging binary pulsars a "universe" 
almost always has its mergers predominantly in spirals 
(dotted curve). 

Similarly, since a significant fraction of observed short 
GRBs are seen at low redshifts (as opposed to during the 
epoch of elliptical galaxy assembly), spiral-dominated 
models usually fit observations better than elliptical- 
dominated ones, as the latter could fit only given an 
unusually long preferred delay between binary birth and 
merger. For this reason, demonstrated by the dashed 
curve in the bottom panel of Figure [Til the "universes" 
which best fit the short GRB redshift distribution - that 
is, which have short GRBs occurring relatively recently, 
during the epoch of spiral star formation - largely pro- 
duce their short bursts in spiral galaxies. 
Low spiral fraction implies long elliptical delays: Con- 
versely, in order for elliptical galaxies to host a signif- 
icant fraction of mergers, those elliptical galaxies must 
produce binaries with an unusually long characteristic 
delay between birth and merger. For example, in the 



right two panels of Figure [T2] compare the prior (top) and 
constrained (bottom) distribution of spiral fraction and 
median delay time i(50%) between birth and merger. As 
seen in the top panel, the most common delay between 
birth and merger in a randomly chosen elliptical popu- 
lation synthesis simulation is around 8 Gyr; however, as 
seen in the bottom panel, for the handful of "model uni- 
verses" which have fewer than 80% of double neutron star 
mergers born in spiral galaxies, the median characteristic 
delay time is at least an order of magnitude larger. 

5.2. Comparisons II: BH-NS as burst source 

On the one hand, for purely technical reasons described 
in §[4] we cannot require our BH-NS "model universes" to 
be consistent with the present-day NS-NS merger rate. 
On the other hand, unlike the NS-NS few (84 out 

of 500) BH-NS "model universes" do predict too few bi- 
nary mergers to reproduce observations, even given the 
most optimistic assumptions regarding beaming and the 
minimum luminosity of short GRBs; see Figure [5J Fur- 
ther, as shown in Figure [UJ a significant fraction of red- 
shift distributions (183 out of 500, or ~ 37%) are con- 
sistent with the limited sample, and many of those (137 
out of 500, or 27%) are consistent with both the short 
GRB redshift distribution and detection rate; see Figure 
[TOl From these constraints, we can draw the following 
conclusions: 

Mergers usually in spirals: Just as with NS-NS models, 
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Fig. 11. — Distribution of the fraction of BH-NS mergers (top 
panel) and NS-NS mergers (bottom panel) expected to occur in 
spiral galaxies at the present day [Eq. II11H evaluated at z = 0]. 
Solid curves include all models; dashed lines include only those 
models reproducing the redshift distribution; dotted line includes 
only those models reproducing the merger rate of NS-NS binaries 
in spiral galaxies, based on the Milky Way; and the dot-dashed 
line includes only models which could possibly produce as many 
short GRB events as are observed. To better indicate a very high 
or very low fraction f s of mergers occurring in spirals, we plot the 
distributions above versus X = tanh _1 (2/ s - 1). 
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Fig. 12. — A mixed contour and scatter plot of the spiral fraction 
f s versus the characteristic delay t(50%) between birth and merger 
in elliptical galaxies, for BH-NS models (left) and NS-NS models 
(right). The top panels include all 500 models; the bottom panels 
provide only those models which satisfy both constraints (for BH- 
NS simulations, the 137 "model universes" which are consistent 
with the short GRBs redshift distribution and detection rate; for 
NS-NS simulations, the 35 "model universes" which are consistent 
with the short GRB redshift distribution and the number of merg- 
ing binary neutron stars in the Milky Way). This figure illustrates 

(i) that the vast majority of constraint-satis fyin g models have most 
mergers in spirals [also illustrated in Figure [Tl] and more critically 

(ii) that the few models which permit most mergers to occur in 
elliptical galaxies are slightly biased towards longer characteristic 
delay times t(50%). 



though a priori our "model universes" are equally likely 
to produce mergers and short GRBs in elliptical and spi- 
ral galaxies, the set of "model universes" which repro- 
duce a short GRB redshift distribution that is dominated 
by recent mergers is biased towards spiral- dominated 
models: compare the solid and dashed curves in the top 
panel of Figure [TlJ However, unlike the NS-NS case 
where both the redshift distribution (preferring recent 
star formation) and double neutron star merger rate (re- 
quiring very high merger rates in spiral galaxies) limit us 
to "model universes" that produce extremely many merg- 
ing binaries, for BH-NS "model universes" only the red- 
shift distribution carries any information that can bias 
results in favor of spiral galaxies [Fig. Ill] , In particu- 
lar, as illustrated by Figure [71 the distribution of spiral- 
galaxy BH-NS merger rates remains nearly the same no 
matter what constraints have been applied. 
Ellipticals can be hosts only with significant characteristic 
delays: Further, as one expects on physical grounds, el- 
liptical galaxies can contain a significant fraction of short 
GRBs and mergers only when unusually long character- 
istic delays between birth and merger are involved; see 
the left two panels of Figure [T2l 

Unlike the double neutron star case, a significant frac- 
tion of "model universes" predict low spiral fractions 
f s < 0.8. In order for old elliptical galaxies to host a 



number of mergers similar to the young spiral galax- 
ies, the binaries in old elliptical galaxies must survive 
over rather long times - many Gyr after their forma- 
tion. Comparatively speaking, population synthesis sim- 
ulations of BH-NS binaries produce more systems with 
the required characteristic long delay time than do sim- 
ulations of NS-NS binaries; see for example the fraction 
of simulations in the top right panels in either Fig. [5] or 
the top left panel of with i(50%) greater than 10 Gyr. 

5.3. What fraction of mergers produce short GRBs? 

So far, we have required our "model universes" to pro- 
duce when given the most favorable assumptions regard- 
ing burst luminosity distribution and beaming at least 
as frequent short GRB detections as are observed, as 
shown in Figure [5J This fairly weak constraint is im- 
posed because, under perfectly plausible but less opti- 
mistic assumptions, a great many short GRBs could be 
missed. However, the difference in Figure [5] between, on 
the one hand, the vertical lines indicating the observed 
short GRB detection rate and, on the other hand, an 
optimistic prediction fdfbRD can be reinterpreted as the 
fraction of short GRBs that must be missed in order for 
our predictions to correspond with observations. 

By way of example, if a NS-NS "model universe" corre- 
sponds to an optimistic detection rate of order 10 yr _1 , 
then only 1 of every 100 NS-NS mergers could produce 
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short GRBs brighter than the least luminous seen and 
aimed in our direction. This factor of 1/100 could be 
produced by any one or combination of several factors, 
all of which act to decrease the predicted short GRB de- 
tection rate below these optimistic predictions: (i) beam- 
ing, since the predicted detection rate is proportional to 
fb and our optimistic calculations assume /& = 1; (ii) 
fainter short GRBs, since the detection rate is also pro- 
portional to the photon luminosity N m { n of the faintest 
short burst (i.e., oc N m - ln / N m i nsccu ); (hi) some intrinsic 
physics which prevents all but a select few mergers to 
produce detectable bursts; or even (iv) further changes in 
our population synthesis model, such as assuming fewer 
than 100% (our present choice) of all stars are born in 
binaries. 

No incontrovertible evidence exists that requires any 
of these factors be less than unity. However, there is 
good observational and theoretical motivation for re- 
evaluating our detection rate constraint using a prior 
on the product of all these factors: essentially, while all 
could be nearly unity, good reasons exist to imagine that 
se veral may be significant ly less than 1. For example, 
(i) iSoderberg et alj (|2006aD and others have argued that 
short GRBs may be beamed. Further, (ii) the least lu- 
minous burst seen is nearly at our detector's sensitivity 
limit, yet short GRBs peak flux distribution is a fea- 
tureless power law. Since only a remarkable coincidence 
could produce such a fortuitous combination of detector 
design and source strength that this burst is indeed at or 
near the intrinsic burst sensitivity limit, we can reason- 
ably expect the minimum-luminosity burst to be signifi- 
cantly less luminous than t he faintest burst yet seen . Ad- 
ditionally, as described in Bclczvns ki et al.l (|2007i), (iii) 
theoretical simulations of BH-NS mergers could possibly 
produce short GRB merger events only for a limited ar- 
ray of binary component masses and spins; the fraction 
of mergers which could be short GRBs ranges from 1/3 
to 1/50, depending on the BH birth spin. Finally, (iv) 
based on comparison with the present-day Milky Way, 
the binary fraction could be u p to a factor two lower 
than the value we assume; e.g., IBelczvnski et al.l (|2007f ) 
adopts a binary fraction of 50%. If each of these four 
factors is decreases the fraction of short GRB detections 
below our predictions by only factor of roughly 3, then 
our expectations for the short GRB detection rate cor- 
responding to a given "model universe" should be lower 
by nearly a factor of 100! 

Adopting a canonical factor of 1/100 to transform the 
optimistic detection rates presented in Figure [8] into a 
"realistic" expectations leads to a dramatic transforma- 
tion of our understanding. On the one hand, extremely 
few BH-NS "model universes" would produce short GRB 
events as frequently as are observed. On the other, these 
"realistic" assumptions cause the 35 constraint-satisfying 
NS-NS "model universes" to predict roughly as many 
short GRB events as are observed (i.e., imagine shifting 
the dotted peak in Figure[5]to the left by 2 orders of mag- 
nitude). While this tantalizing correspondance could be 
a coincidence and while the precise results of our "realis- 
tic" prior cannot be taken too seriously, they do remind 
us of two salient features of our predictions: (i) that BH- 
NS "model universes" are consistent with the data only 
given relatively optimistic assumptions, and that they 



quickly become inconsistent with observations if those 
assumptions are relaxed; and (ii) that on the other hand 
because NS-NS "model universes" can allow for so many 
detections and therefore have much more flexibility in 
the fraction of short GRBs that are missed, and because 
double neutron star models are otherwise entirely consis- 
tent with all existing observations, double neutron star 
models for short GRBs remain an attractive candidate 
explanation for short GRBs. 

Comparison to related studies: IBelczvnski et all (|2007t ) 
(hereinafter BTRS) also used the StarTrack popula- 
tion synthesis code and hydro dynamical simulations of 
BH-NS mergers with BH spins ((Rantsiou et al.ll2007| ) to 
determine whether BH-NS merger events occurred fre- 
quently enough, with the right combinations of parame- 
ters, to reproduce short GRB observations. In contrast 
to the present broad study, which relics only on observa- 
tional constraints, this study adopts several of the "re- 
alistic" priors mentioned above to reduce the fraction 
of merger events that produce short GRBs. Specifically, 
this paper relies on well-motivated hydrodynamical stud- 
ies to argue that at best a small fraction of BH-NS merg- 
ers (1/3 to 1/50, depending on black hole birth spin) will 
produce burst events. Additionally, the BTRS simula- 
tions rely on a single "most-plausible" set of population 
synthesis parameters, including in particular a high com- 
mon envelope efficiency aA 15 . As a result, compared to 
the wide range of common-envelope efficiencies used in 
this study, BTRS find relatively low BH-NS merger rates. 
We also note that BTRS assume a 50% binary fraction, 
lower than the 100% used here. Combining these three 
factors, 16 BTRS conclude that BH-NS mergers occur too 
infrequently (less than 10Gpc _3 yr _1 ) to explain short 
GRB merger rates. 

6. CONCLUSIONS 

In this paper, we used a large archive of concrete, cur- 
rent population synthesis calculations to generate merger 
rate densities for NS-NS and BH-NS binaries. Using as- 
sumptions regarding short GRB source luminosities and 
detection selection effects, we then compared these rate 
densities to short GRB detection rates and redshift dis- 
tributions, as well as to (when available) the present-day 
Milky Way binary neutron star merger rate. Whether 
assuming short GRBs arose from either BH-NS mergers 
or NS-NS mergers, a small but still substantial fraction 
of models were consistent with existing observations us- 
ing the most optimistic assumptions for certain priors 
(i.e., no beaming and all s tars born in binari es). Con- 
trary to an earlier study bv lNakar et alj (|2005l ). we have 
demonstrated by using a two- component star formation 
model that exceptionally long characteristic delays be- 
tween binary birth and merger are not uniquely required 
needed to reproduce the short GRB redshift distribution 
and detection rate, though they are strongly preferred 
if we additionally demand a significant fraction of short 

15 The common envelope efficiency reflects the fraction of orbital 
energy needed to eject the envelope; since most BH-NS binaries go 
through a common-envelope phase, a low efficiency implies these 
binaries will have a tight final orbit. 

16 Additionally, the simulations used in this paper allowed much 
less mass accretion onto black holes during the common-envelope 
phase. Based on studies conducted in our group, this change does 
not produce a dramatic difference in merger rates. 
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GRBs occur in elliptical hosts. Generally our simula- 
tions favor spiral- dominated models: f s > 0.7. Double 
neutron star models reproduce observations only if spi- 
ral hosts dominate; BH-NS models can permit a wide 
range of spiral fractions. Though we have not imposed 
it as a constraint, the fraction of short GRBs in spiral 
hosts potentially provides an extremely powerful addi- 
tional constraint on source models. For example, the 
six host identifications shown in Table [1] suggest the ob- 
served spiral fraction f s is 50%. However, to impose it 
reliably requires a careful study of the systematic uncer- 
tainties in the two-component star formation model used 
as well as of the selection effects in host galaxy identifica- 
tions. For example, more massive elliptical galaxies will 
more efficiently retain any strongly-kicked neutron star 
binaries but have less gas into which a the short GRB 
blast wave can collide and produce an afterglow; these 
and other selection effects 

Generally speaking, the exceedingly small number of 
well-identified short GRBs strongly limits our ability to 
draw conclusions based on comparisons to properties of 
that set, such as to the redshift distribution or to the frac- 
tion of bursts seen in spiral hosts. More well-identified 
hosts are needed; additionally these hosts should hope- 
fully be drawn from a less -biased sample than the present 
sample appears to be (see lBerger et al-lfeOOGl ). Of course, 
we encourage deep follow-up searches for afterglows, par- 
ticularly since our model universes almost always predict 
a higher proportion of high-redshift short GRBs than 
has yet been observed. We particularly encourage the 
development of thorough redshift-limited surveys, since 
a determination of the relative proportion of bursts in 
star-forming or elliptical hosts can be done in the nearby 
universe has significant potential to improve our under- 
standing of the formation mechanism. 

Though the limited sample of well-characterized short 
GRBs currently limits our ability to constrain input 
physics, we expect that the significant increase in 
event statistics expected over the next few years will 
make these events a primary mechanism to constrain 
double compact object merger rates and the associ- 
ated astrophysics. For example, with well-understood 
short GRBs already outnumbering the few known dou- 



ble neutron stars in our galaxy, short GRBs would 
soon provide the most precise (nonparametric) obser- 
vational constra int on compact object merger rates (cf 
iKim et afll2003f ). provided the selection effects are quan- 
titatively understood. In turn, th ese constraints inform 
us about binary stellar evolution (jO'Shaughnessv et alJ 
(2007bT) [Beiczynski et al.l (|2006ah ). Additionally, since 
short GRB sources would also be copious emit- 
ters of gravitational waves, Swift and ground-based 
gravitational- wave observatories operating in coincidence 
could conceivably probe th e details of an unusually close 
merger event itse l f (see,e.g.,[k obavashi & Mes zarosl2003t 
iFinn et all [20041 : iNakar et all l2005h. LIGO in particu- 
lar is both already operating at design sensitivity and 
conducting triggered searches from GRB observations 
([Abbott fc et alH2005h . Given its compelling qualitative 
agreement and far-reaching scientific impact, the merger 
hypothesis deserves a detailed comparison against the 
best possible models, to either corroborate it, definitively 
disprove it, or discover weaknesses in our assumptions 
and models. 

Finally, the current sample does not yet allow us 
to quantify the minimum luminosity and beaming of 
short GRBs. However assumptions regarding the frac- 
tion of merger events which do not produce short GRBs 
brighter than the weakest burst seen so far and point- 
ing towards us combin ed with theoretical priors (see, e.g. 
Bclc zvnski et al.| [2007). indicate that BH-NS mergers do 
not occur frequently enough to explain all short GRBs. 
We anticipate that future larger samples will allow us to 
place stronger constraints on the merger models. 
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APPENDIX 
ARCHIVE SELECTION 

The predicted short GRB detection rate depends directly on the fraction of star forming mass A predicted to end 
up in GRB progenitors. However, the mass efficiency A can vary substantially depending on model assumptions. In 
order to constrain the range of short GRB detection rates expected (for a fixed minimum short GRB luminosity N) , 
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Fig. A13. — For spiral (left two panels) and elliptical (right to panels) archives, a scatter plot of the number of NS-NS (top) and BH-NS 
(bottom) binaries n seen in a simulation of N progenitor binaries (all of m > 4A/q , with primary mass drawn from the broken-Kroupa 
IMF). The two solid lines show the cutoffs n > 20 and nN > 5 X 10 6 imposed to insure data quality and reduce sampling bias. Simulations 
used in this paper, shown in black, must lie above and to the right of these cutoffs. 



we require the unbiased distributrion of A. However, the large variety of simulation sizes in our simulations introduces 
a bias in our estimate of the mass efficiency distribution - or, equivalently, in the distribution of n/N for n the number 
of binaries seen and N the number of binaries simulated. Specifically, not all of our archived population synthesis 
simulations contain enough of each type of event (indexed by K) to produce reliable predictions involving it. The set of 
simulations with n greater than any fixed threshold threshold (including n — 0) is biased, over-representing simulations 
with high n/N; see for example Figure IA13I Since the mass efficiency A is directly proportional to n/N [Eq. ([6])], we 
describe a filter on n and N which preserves the distribution of n/N (for the distribution of A) and simultaneously 
insures that the average simulation has many binaries ((n) > 200, so each model's delay time distribution dP/dt can 
be accurately estimated). 

Formally, each population synthesis simulation converts some unknown fraction s of progenitor binaries into target 
events n — sN. In practice our population sample size N is roughly randomly chosen, independent of the unknown s. 
In other words, we expect the distribution of n and N to derive from two independent distributions for N and s, with 
densities / and g: 

dP(\nN,\nn) = f(\nN)g(\n(n/N))d\nNd\nn 

= f(\nN)g{\nn/N)d\ny/Nnd\nn/N (Al) 
The distribution of s — n/N can be quickly extracted from any two-dimensional distribution in n,N by: 

gdlns = j ' 8{\n{s/{n/N)))dP . (A2) 

However, the above method assumes the distribution perfectly resolved. In practice simulations are generally not 
repeated, so "fractional" and small n cannot be resolved by repeated trials; instead, only simulations with s > 1/f 
will produce enough events to provide reliable estimates of s, and thus its distribution. This truncation effect biases 
our reconstruction of the two-dimensional density and therefore our reconstruction of the mass efficiency distribution. 
If we had perfect foreknowledge, however, we could have chosen our sample size so nN was constant and large. Though 
we would choose n and TV in a highly correlated fashion, we would guarantee each point was well-resolved; our method 
above would correctly reconstruct the distribution of s. Hence choosing all population synthesis archives with nN 
greater than any threshold will allow accurate estimation of the distribution of s. 

Based on the above discussion and the observed distributions in n and N, we introduce cutoffs on n (> 20) and 
nN (> 5 x 10 6 ) which remove the least resolved simulations from consideration (via the n cutoff), reduce the bias 
associated with a minimum n (via the nN cutoff) , and insure that most simulations have at least 200 binaries (both) . 
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Fig. B14. — Results for errors in estimates of P(t m ) and dP m /dt, versus the number of points used in the estimate, for the trial 
distribution function dP m /dt = °'^ 79 between 30 and 10 7 Myr. Left: Plot of the average and 90% supremum-norm error versus n for 
Pm(t), along with fits to these two quantities (namely, n _0,45 /3 and 0.7n~ ' 47 , respectively). Right: Plot of the average and 90% bound on 
fractional error in dP m /dt versus the number of points available, along with fits to these two quantities (namely, 0.36n -017 and 0.73n — 023 , 
respectively). 



ESTIMATING DELAY TIME DISTRIBUTION 

Poisson errors associated with the limited sample size of our population synthesis simulations inhibits our ability 
to reconstruct the delay time distribution, whether represented as a cumulative distribution P m {t) = P m (< t), the 
probabili ty that a delay between birth and merger is less than t, or as dP m /dt. Using classical statistical methods - 
see, e.g., iMerrittl (|1994f ) . iThompson fc Tapial (|1976( ), and references therein - we smooth the n observed delay times 
to build our estimates for P m (t) and dP m /dt. This appendix briefly reviews those methods and the accuracy of the 
resulting estimates. 

Estimating cumulative distribution: We estimate the cumulative distribution P(< i) for a sample of population 
synthesis events by a function P(t) that smoothes sample dataover a short smoothing length s in l t = log 10 (£/Myr): 

n 

P s (l t ) = J2®s(k-k,k) (Bl) 

fe=l 

e s (x) = - (l + eif(x/sV2)) (B2) 
[(max fc Z tjfc ) - (min fc Z tjfe )] 

Sl = io^ (B3) 

for n the number of events in our population synthesis sample, k = 1 . . . n an index over tat same sample, and 
It.k = log(i/c/Myr) the logarithm of the delay tk between binary formation and merger for the fcth binary. The 
extremely short smoothing length used here approximately minimizes the average difference between predictions and 
results. 

To test this approach, for several n we drew many Monte Carlo samples of n delay times from a canonical cumulative 
distribution: 

Pm{k) = lt 7 ZtS Q{lt ~ 1O g( 3 °)) ( 7 " l ^ 

Figure IB 141 shows the average and 90% distance max/JP m (/ t ) — P m (lt)\ between our function and our fit versus the 
number of points smoothed to estimate P m (t). 

Given the archive selection procedure presented in ^ each archive typically contains of order n ~ 100 merging 
binaries. For such a typical archive, our smoothing method will reconstruct P(t) almost everywhere to better than 5%, 
with the largest errors typically arising near the largest and smallest delay times. As expected, the maximum error 
agrees perfectly with the distance between observations and data for a Kolmogorov-Smirnov 95% hypothesis test. 

Estimating differential distribution : Similarly, we estimate dP m /dt by P mi io s , defined using the time derivative of 
Pm{< t) after converting from logarithmic to physical time: 

p = 1 dPm - 1 V 1 p -(\ogt-l t , k ) 2 /2(s 2 f (m) 

m UnlO dl t tin 10^ V27r(s 2 ) 2 V ' 

[(maxfc/f.fc) - (mm k lt,k)] , oc x 

s 2 = : - 7= — (B5 

1.25-v n 
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Compared to the previous case of estimating the cumulative distribution P m (t), where a significant fraction of all 
sampled points affect the estimate at any t, a significantly longer smoothing length S2 is required to estimate dP m /dt, 
since roughly only those few sample points within S2 of t contribute to our estimate. As a result, we cannot confidently 
estimate dP m /dt more accurately than w 30%. 

To demonstrate that smoothing produces a fairly inaccurate prediction for dP m /dt, we applied it to the trial problem 
mentioned above. To estimate the mean relative error / associated with our estimate in the physically pertinent interval 
of roughly 100 Myr to 15 Gyr by 

1/2 

(B6) 



4.5 



dl t \P-P\'- 
2 - 7 (P 2 + P 2 )/2_ 

which generalizes an rms measurement of the relative error in P to allow for large relative errors [i.e., when P m — 
P m (l + S(l t )), then I = (J 5 2 dl t / Alt) 1 / 2 ]. The second panel of Figure [BM] demonstrates that large errors are inevitable, 
even with many sample points and using a smoothing length S2 which is appr oximately optimized for ea ch n. While 
this slow convergence is a familiar problem for all density estimators (see, e.g. iThompson fc Tapial 1 1976t h it severely 
limits the ability of any population synthesis simulation to reliable estimate dP m /dt, given the severe computational 
limits involved. By way of example, a delay time distribution accurate to 5% would require roughly 6 4 ~ 1300 times 
more binaries than a 30% accurate estimate. 

Implications for short GRB predictions: The short GRB detection rate, redshift distribution, and elliptical-to-spiral 
ratio all depend on estimates of dP m /dt. However, while our reconstruction of dP m /dt for any fixed population 
synthesis model involves substantial uncertainties, these are not much greater than our corresponding uncertainty in 
our estimate of that model's the mass efficiency A. Further, these uncertainties are vastly less than the systematic 
uncertainty involved in not knowing the physically appropriate population synthesis model (i.e., roughly two orders of 
magnitude uncertainty in A). 



